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ABSTRACT 

We investigate the orbital structure of a model triaxial star cluster, centered 
around a supermassive black hole (BH), appropriate to galactic nuclei. Sridhar and 
Touma (1999) proved that the presence of the BH enforces some regularity in the 
dynamics within the radius of influence of the BH. We employ their averaging method 
to reduce the degrees of freedom from three to two. Numerical orbit integrations, 
together with Poincare surfaces of section allow us to draw a global portrait of the 
orbital structure; in our calculations we employ a model cluster potential that is 
triaxial and harmonic. The averaged dynamics of the axisymmetric case is integrable, 
and we present a detailed comparison of orbits in oblate and prolate axisymmetric 
potentials. Both cases support resonant orbits with fixed values of eccentricity, 
inclination, and periapse, whose lines of nodes rotates steadily. We then systematically 
explore significantly triaxial potentials, possessing small oblateness, or prolateness. 
Resonant orbits and their families are studied both numerically, and through secular 
perturbation theory. Chaos appears to be suppressed for all the cases we studied, and 
we obtain effective third integrals. Some of the orbits appear to reinforce the shape of 
the potential; we provide phase space, as well as real space portraits of these orbits. A 
particularly promising resonant orbit exists in highly prolate, triaxial potentials. 



Subject headings: black hole physics - celestial mechanics, stellar dynamics - galaxies: 
nuclei 
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1. Introduction 

Evidence for the presence of central black holes (BH) of mass, 10^ — 10^'^ Mq, appears strong 
for about a dozen galaxies (c.f. Kormendy and Richstone 1995). The masses of these objects are 
consistent with the mass in BHs needed to produce the observed energy in light from quasars 
(Soltan 1982, Chokshi and Turner 1992). Recent studies suggest that the dynamical influence of 
the BH extends to well outside the nuclear regions of its host galaxy, if the orbits of stars carry 
them close to the center (c.f. Merritt 1998 for a review). In particular, dynamical phenomena in 
the central regions are strongly influenced by the BH. Stars move in the combined gravitational 
potentials of the BH, and the surrounding mass of stars and gas (for brevity, these will be referred 
to as the "cluster"). The radius of influence of the BH (r/j) may be defined as radius of the sphere 
within which the mass of the surrounding stars and gas is less that of the BH. Within r^, the 
potential of the BH dominates the dynamics, and the potential of the cluster may be considered to 
be a small perturbation. Sridhar and Touma (1999; hereafter referred to as ST99) studied stellar 
dynamics within rh, and arrived at the following conclusions: 

(i) Stellar orbits may be thought of as Keplerian ellipses that deform and precess over times that 
are longer than orbital times by a factor (M/Mg), where Mg is the mass of the cluster enclosed by 
the orbit. Well within rh, M/Mc » 1. 

(ii) The slow dynamics of precessional motions may be understood by averaging over the Keplerian 
orbital periods of stars, which are the shortest time scales in the problem. 

(iii) Averaging gives rise to a secular integral of motion, / ~ y/GMa, where a is the semi-major 
axis of the orbit. Thus the BH enforces some degree of regularity in the structure of orbits that 
lie within rh- 

(iv) A stellar orbit in a time-independent cluster conserves two integrals of motion. These are the 
energy, which is an exact invariant, and the semi-major axis, which is secularly conserved. 

(v) Thus clusters which may be modelled as razor-thin disks, or as three dimensional, axisymmetric 
objects, give rise to averaged dynamics that is integrable. 

(vi) Non-axisymmetric, razor-thin disks, with power law surface density profiles have averaged 
dynamics, whose qualitative features are independent of the steepness of the power law. This 
occurs because, within rh, the Kepler potential of the BH is far steeper than any self-gravitational 
force due to the cluster. 

Here we take up the thread, and consider orbits in three dimensional clusters, which is a far 
more complicated problem than the cases studied by ST99. We specifically wish to focus on the 
effects of three dimensionality, and triaxiality. Recalling from item (vi) above that the qualitative 
features of orbital structure (arrangement of lens and loop orbits) are insensitive to the steepness 
of the cusp, we wish to choose a cusp steepness that would yield the simplest expression for the 
averaged Hamiltonian for our three dimensional problem. For the razor-thin disks studied by 
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ST99, the averaged Hamiltonian turned out to be a hypergeometric function of the eccentricity, for 
a general value of the cusp steepness; on the other hand, the harmonic case yielded the simplest 
expressions, composed only of elementary functions. Therefore, in our present venture into three 
dimensions, we restrict attention to a family of triaxial, harmonic potentials, average over the 
fast orbital motion, to obtain an expression for the Hamiltonian governing slow dynamics. As 
equations (4) and (5) show, only elementary functions occur in this Hamiltonian, which governs 
nearly Keplerian, three dimensional dynamics. The phase space of the averaged dynamics is four 
dimensional, hence the dynamics is well suited to exploration through Poincare surfaces of section. 
Because of all these favourable features we consider the Hamiltonian of equations (4) and (5) as 
the "canonical" Hamiltonian governing three dimensional dynamics within the radius of influence. 
However, these nice features may be invalid when the orbital radius becomes comparable to r/j; 
in particular, the discrepancy between true dynamics, and averaged dynamics may be greater in 
three dimensions. 

In § 2, we provide a quick introduction to action-angle variables for the Kepler problem, 
and averaging over the fast orbital frequency. We then study the orbital structure as a function 
of the energy (after subtracting out the Keplerian contribution). Averaging over the Keplerian 
motion, promotes the semi-major axis as a conserved quantity. For scale-free potentials, such 
as the harmonic potential, the orbital structure is similar at all radii for which the averaging 
procedure is valid; this motivates the scalings adopted later. Orbits that reinforces the shape of 
the generating potential, demands special attention. Such orbital families play a crucial role in the 
construction of self-consistent stellar dynamical models (e.g. Schwarzschild 1979, de Zeeuw 1996, 
Merritt 1999). Thus we also study the real space structure of orbits, when they appear worthy of 
mention. Averaged dynamics for oblate and prolate axisymmetric cases is discussed in § 3. We 
numerically integrate the equations of motion, and take Poincare sections. In § 4, we discuss the 
global orbital structure for the triaxial case, when departures from spherical symmetry are small 
(small oblateness or prolateness, but quite large triaxiality). The results are interpreted through 
secular perturbation theory. 

The reader who is interested in getting a quick glimpse of the real-space structure of those 
orbits that reinforce triaxiality, as well as oblateness/prolateness of the potential, is advised to 
consult the list of relevant figures and comments provided at the end of the discussion in § 5. 

2. Formulation of orbit— averaged dynamics 

Let us locate the BH, of mass M, at the origin (this is a valid assumption, so along as the center 
of mass of the surrounding cluster does not vary with time). We introduce action-angle variables, 
appropriate to the Kepler problem, that allow us to conveniently carry out orbit-averaging. These 
are the Delaunay variables (c.f. Plummer 1960, Goldstein 1980). / = ^/GMa, L, and are 
the actions, where L and are the magnitude and z-component, respectively, of the angular 
momentum. We denote the conjugate angles by w, g and h, respectively, w is the orbital phase 
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Fig. 1. — Orbital elements of the Kepler problem (Delaunay variables): the orbit is represented by 
the ellipse, drawn in bold line. P marks the pericenter. The figure shows only the four variables of 
precessional dynamics. The other two refer to the semi-major axis of the ellipse, and orbital phase. 
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("mean anomaly") which varies from at pericenter to 2tt after one orbit, g is the angle to 
the pericenter, measured from the ascending node, and h is the longitude of the ascending node 
(see Figure 1). For the Kepler problem, the only Delaunay variable that is time dependent is w: 
w{t) =wo + {G^M^/I^)t. 

Let ^{x,y,z) be the perturbing gravitational potential of the cluster. To average $ over the 
orbital phase, it is convenient to express the spatial coordinates in terms of the Delaunay variables 
(c.f. Goldstein 1985): 



/ X \ 



/ CqCfi — CiSliSq 



V ^ / 
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V 



SiSn 
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(1) 



where S and C are shorthand for sine and cosine of the angles given as subscript, i is the angle of 
inclination determined by cosi = Lz/L, and e = y^l — L'-^/P is the eccentricity; as before, a is the 
semi-major axis, r] is the eccentric anomaly, related to the orbital phase through w = r] — es'mr] . 
Having substituted the expressions for {x,y,z), given in equation (|l]), in ^{x,y,z), averaging is 
conveniently carried out by integrating over ij: 



27r 



dr]{l — e cos r/) $ . 



(2) 



This orbit-averaged potential, ^, plays the role of the Hamiltonian for slow, precessional dynamics 
(see ST99). Superimposed on the slow time variation of averaged dynamics are fast oscillations 
(with frequency of order ^/GMja^), whose fractional amplitudes are of order e ~ a ($ — ^)/GM . 
Henceforth we assume that e <C 1, and ignore these fast oscillations. Note that, for fixed a, e is 
larger for more eccentric orbits, because these sample a larger spatial variation of $ . Because <I> is 
independent of w, the conjugate variable, /, is secularly (i.e. on the average) conserved. 

For investigation of orbital structure, we consider a family of triaxial, perturbing potentials 
of the form 



<^>{x,y,z) 



+ 



(3) 



'0 \ ^ip/ 

Physically, this gravitational potential may be assumed to arise from a homogeneous, triaxial 
cluster whose center coincides with the location of the BH. The shape parameters of the potential 
(viz. and c^) are related to those of the density (say, bg and Cg) and can be found in the 
literature (c.f. Chandrasekhar 1969, Binney and Tremaine 1987). The case when the potential is 
axisymmetric may be studied by setting b^p = 1; oblateness, or prolateness is achieved by choosing 

to be less than, or greater than unity, respectively. More realistic, inhomogeneous mass 
distributions could be considered, but at the expense of considerably greater algebraic complexity. 

Let /i = {^qu'^ / Ivq) be a measure of the precession frequency for orbits with semi-major axis 
a. As in ST99 , it proves to be convenient to define a dimensionless time r = fit, Hamiltonian 
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H = (<I>/^I), and angular momenta i = L/I and £z = Lz/I. Note that < £ < 1, and 
— 1 < ^2 ^ 1 • Averaging the potential of equation (|3|) over orbital phase, we obtain the 
dimensionless Hamiltonian 

Hii, 4, g, /i) = 1(5 - 3f) + eM£, 4,5) + eM£, 4, 9, h) , (4) 

where 

Hc{£,£.,g) = \(^l-^^^{5-3£^-5{l-£^)C2g) , 

Hb{£,£z,g,h) = (^^-2£^^(^CgSh + jChSgy + ^^(^-SgSH + jChCg^\ (5) 

and Efe = (6^^ — 1), Ec = (c^^ — 1). This is the model two degree-of-freedom Hamiltonian that will 
be studied in the rest of this paper. The equations of motion are, 

■ dH dH ■ dH ■ dH 

where the over-dots refer to {d/dr) . We integrate these equations numerically using a fourth- 
order, adaptive step-size, Runge-Kutta scheme, and take surfaces of section at constant values 
of either of the angles g, or h. Thus we follow the the deformations of Keplerian ellipses of fixed 
semi-major axes, over time scales much longer than the orbital periods. This technique will, 
in principle, reveal all of the phase space structure of the orbits. We will also supplement the 
numerical calculations with analytical estimates, wherever necessary. Note that such an assault on 
orbits in triaxial potentials would have been impossible, without the presence of the extra integral 
of motion that averaging promotes. 

The Hamiltonian, given in equations (^) and (^, is a constant of motion. The orbital 
structure shows strong dependence on the value of this quantity, and henceforth we will refer to it 
as the energy (E). Since, for the averaged dynamics, Kepler energy is constant, this refers refers to 
the average value of ^ over the orbit in question. The simplest case is that of spherical symmetry, 
when eb = ec = 0. Then H = {5 — 3^^)/2, so that £ = £z = h = , and g = —3£ : the Keplerian 
ellipses precess rigidly, at a constant rate in a retrograde sense, maintaining an invariant orbital 
plane. 



3. Orbits in the axisymmetric case 

When bip = 1, the potential of equation is axisymmetric about the z-axis. Hence is 
conserved even for the unaveraged dynamics. For slow dynamics, this implies that £z is an integral 
of motion, and we are left with a one degree-of-freedom system, whose dynamics is obviously 
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integrable. Setting eb = 0, from equations @ and (|5|), we obtain the following expression for the 
axisymmetric Hamiltonian, 



(7) 



where iz is now regarded as a constant parameter. Note that €c is positive for oblate configurations, 
and negative for the prolate ones. The nodal precession rate given by 



shows that the sign of h is opposite to that of the product, ec^z; hence nodal precession is prograde 
(retrograde) for prolate (oblate) potentials. 

To obtain a global picture of the dynamics, it is instructive to plot an £ — g (Poincare) surface 
of section taken at some fixed value of h; Figures 2a and 2b show two such sections for oblate 
(c<^ < 1) and prolate (c,^ > 1) potentials, respectively. 

From these figures, it may be seen that the dynamics is organised around the elliptic fixed 
points, located at g = tt/2 and 3it/2. These are special orbits for which i and g are constant in 
time, whereas h precesses steadily. Using the Hamiltonian of equation (0) in equations (^), it is 
straightforward to verify that these fixed point orbits (FPOs) are described by. 



g = 7T/2 or 3^/2, h = ieJz-5eJ Z ■ (9) 



Before we proceed to examine the dynamics in more detail, let us pause to consider some 
implications of equations (^) for FPOs (it is useful to recall that 0<£<1,— 1<£2<1, and 
cos i = iz/^)- 

1. One end of the FPO sequence is set by £z = 0, which is the z axial orbit (unless 3 + 4ec = 0, 
which is the case c^, = 2). 

2. Because oc cos^ i, the orbits become rounder as their inclination decreases (i.e. as the orbital 
plane approaches the equator). 

3. The other end of the FPO sequence is reached when either i = 0, or ^ = 1, whichever occurs 
first. This depends on whether 5ec/(3 + 4ec) is less than, or greater than unity (equivalently, when 

is greater than, or less than 1/2). There is an additional requirement that 5ec/(3 + 4ec) > 0, 
which implies that no FPOs of the kind allowed by equation (P) exist, for prolate potentials which 
have 1 < < 2 . 

4. Oblate: < < 1/2: ^ = 1 is reached before i = 0, so there are no equatorial FPOs. 

5. Oblate: 1/2 < < 1: i = is reached before £ = I . When c^p = 1/2, equation (^) implies that 
ell'^ = cos^ i for the FPO; thus both i = and £ = 1 are reached simultaneously. 
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6. Prolate: c^p > 2: £ = 1 is reached before i = 0, so there are no equatorial FPOs. 

The Poincare section may also be regarded from a different perspective. For the completely 
integrable dynamics of the axisymmetric problem, one may invert equation (|^), to express iz as a 
function of i, g and 

We may now regard H as a, constant parameter, with the i — g dynamics controlled by the 
Hamiltonian equal to {—iz), and "time" measured by h. Figures 2a and 2b may be reproduced by 
plotting the isocontours of iz in the i — g plane. It is straightforward to verify the stability of the 
orbits of equation (P) by Taylor expanding iz about the fixed point values of i and g. The FPOs 
parent a family of orbits for which both i and g librate periodically, and h precesses non-uniformly 
in time. In real space, these orbits fill ( "hollowed-out" ) conical regions whose axes coincide with 
the z-axis. The phase flows in the i — g plane are structurally similar to those of a pendulum. 
Hence there are also other orbits for which g circulates, and these are distinguished from the 
former by a separatrix. Figures 2c and 2d show two representative sections of such orbits. 

For prolate potentials with 1 < < 2, the dynamics is quite different, as is revealed by the 
surface of section, shown in Figure 2e; g circulates for all orbits. This difference in behaviour 
is connected with the stability of the z-axis orbit. Let us test the z-axis orbit for stability to 
eccentricity variations which, however, maintain iz = 0. We Taylor expand the right side of 
equation (0) about i = 0, and g = 7r/2. Setting g = it/2 + 6g, we obtain 

^-{l + e,)-H = {^ + 2e,)i' + ^e,{6gf. (11) 

The energy of the z-axis orbit is H = 5(1 + ec)/2, which is the maximum allowed energy in 
the oblate case. Therefore, the left side is always non negative, as are the coefficients of i'^ and 
(Sg)'^; thus the z-axis orbit is always stable. For the prolate case, the energy of the z-axis orbit, 
H = 5(1 + ec)/2, is the minimum allowed energy. Hence the left side is always non positive. The 
coefficient of (Sg)'^ is negative, so the issue of stability depends on the sign of the coefficient of i"^. 
It is clear that the z-axis orbit is stable for €c < —3/4, and unstable when —3/4 < ec < 0. 



4. Orbits in the triaxial case. 

The shape of the potential is given by two axes ratios, and c<^. The definition of the 
Delaunay variables gives a privileged status to the z-axis, so we let this be the axis of symmetry 
for the study of orbits in the axisymmetric case. In particular, we set = 1, and chose c^p < 1 and 

> 1 for oblate and prolate cases respectively. Whereas this is the most convenient choice for 
the dynamics, it implies that the ordering of the lengths of the x, y and z axes are different in the 
oblate and prolate cases. This feature carries over to triaxial configurations as well. Hence we allow 
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and to take any positive value, and define the "triaxiality" to be T = (of — /{of — 03), 
where (01,02,03) are equal to the set (1,6<^,C(^), after the latter has been rearranged (if necessary) 
in descending order of magnitudes. With this convention, for axisymmetric configurations, T = 
for the oblate case, and T = 1 for the prolate case. 

The equations of motion (eqn. ^), resulting from the Hamiltonian of equations (H) and (pf), are 




i 



4 



(1 - mCgSh + -iChSg){-SgSh + ^ChCg) 



+66 



(5 - U''){CgSh + -iChSg)(CgCh " ^ S g) 



^"^{ — SgSh + -jChCg){SgCh + -jShCg) 



.(5 _ Sf _ 5(1 _ f)C2g) + £{1 - p{5C2g " 3) 
4 



{CgSh + —ChSg){4£CgSh + 5j^ChSg) + £SgSh{ — SgSh + —ChCg] 



--e,-±{f>-M^-f>[l-t)C2g) 



(5 _ Ae)-ChSg{CgSh + iChSg) + lChCg{-SgSh + iChCg) 



(12) 



We integrate these numerically, and take Poincare sections by strobing in either of the angles g, 
or h. In this paper, we restrict our study to triaxial configurations that are close to spherical; 
thus h^p and c^p are close to unity (equivalently, e?, and ec are close to zero). This allows us to use 
perturbation theory to understand the results of our numerical computations. 

As in the axisymmetric case we study the orbital structure as a function of the (constant 
value of) Hamiltonian, which we refer to as the "energy" . Recall that the Hamiltonian is equal to 
the (scaled) potential of the cluster. The scaling simply means that we have chosen to explore the 
dynamics of nearly Keplerian orbits with unit semi-major axis. Hence the lowest energy orbits are 
circular, and the highest energy orbits very eccentric. Below we explore both oblate and prolate 
triaxial cases. As representative cases, we consider {h^p,c^p) equal to (0.99,0.96) and (0.99,1.04), 
for the oblate and prolate triaxial potential. It may be verified that the triaxiality is quite large, 
namely T ~ 0.254, and 0.804, respectively. 



4.1. Low energy orbits; i <^ 0{l) 



Recall that iz was a constant of motion in the axisymmetric case; a Poincare section, strobed 
in g would have given horizontal rows of dots in the iz — h plane. Hence this surface of section is 
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very useful in revealing the effects of triaxiality, as may be seen in Figures 3a and 3b, which show 
a pendulum-like dynamics of the iz — h motions. The figures do not, of course, reveal that the 
periapse angle, g, always circulates (but this has been verified by plotting i — g sections, and by 
following the time evolution of g). A glance at the equations of motion ( [l^ ) shows that \g\ ~ 0(1), 
whereas \h\ ~ 0(ec) <C l^]. Therefore, it is useful to imagine the motion of the Keplerian ellipses 
on two well-seperated timescales. At nearly constant values of h and i (i.e. on a nearly fixed 
orbital plane), the precessing ellipse fills a circular annulus. On a longer timescale, the orbital 
plane itself either circulates or librates in h, while oscillating in inclination. As in the axisymmetric 
case, oblate triaxial configurations give rise to retrograde precession of h while prograde precession 
results in the prolate case. In Figure 3a, the librating and circulating orbits are long-axis and 
short-axis tubes, respectively; in Figure 3b, the roles are reversed. These orbits are also found in 
more general triaxial galactic potentials with constant density cores. Their persistence in the case 
when a black hole is present may be understood qualitatively by noting that these tube orbits are 
round {£ ^ 0(1))) so they are not destabilized by the black hole (note that the lower i orbits of 
§ 4.2 are completely different). 

The dynamics can be understood by averaging over the fast angle g. Equations (|^ and (^, 
when averaged over g give 

for the Hamiltonian describing the slower dynamics of £z and h; it is straightforward to verify that 
the isocontours of this Hamiltonian reproduce the surfaces of sections shown in Figures 3a and 3b. 
In this approximation, i emerges as a secular invariant (i.e. no secular variation in eccentricity), 
and the (^-averaged dynamics reduces to a one degree of freedom system in the variables, iz and 
/ij^ The orbits corresponding to circulations of h (short and long axes tubes, for oblate and prolate 
cases, respectively) have larger \£z\, so are more nearly equatorial. The orbit fills a region of space 
that is toroidal and non axisymmetric. The approximate axis of symmetry is along the z-axis, 
and the torus itself could be thought of as somewhat pinched along the intermediate axis {y for 
oblate and x for prolate). 

The orbits for which h librates (and iz flips sign), have smaller values of \iz\, so these are more 
nearly polar (long and short axes tubes, for oblate and prolate cases, respectively). Let us first 
consider the parent polar orbit with iz = 0. For the oblate/prolate cases, the orbit lies in a plane 
perpendicular to the long/short axis. Within the plane, g circulates rapidly and nearly uniformly, 
while i shows equally rapid, small amplitude oscillations; hence the orbit fills an elliptical annular 
region. Figure 4 shows projections of a polar orbit on the three principal planes for a prolate case. 
The orbit lies essentially in the plane of the long and the intermediate axes, so reinforces the shape 
of the potential to some extent. 




(13) 



A more systematic approach to averaging over the fast angle uses first order perturbation theory, whose primary 
advantage over equation (Fsh is that the 0(ec) oscillations in £ are well represented. 
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4.2. High energy orbits: I ~ O(y^) 

High energy orbits are characterised by low angular momentum values, hence are all quite 
eccentric. In general, I ~ 0(-\/|ec|). This regime of energy admits a certain scaling, that allows 
us to parametrise the dynamics by the single quantity; B = eb/|ec|. Let us introduce the scaled 
variables, 

1 5 



/ = ^, lz = ^, K = — [H--], s = ^\ec\T, (14) 
e-A \/\er\ Cc V ^, 



where {l,lz) are conjugate to {g,h), and K and s are the new Hamiltonian and time, respectively, 
for this very slow dynamics of highly eccentric motions. Substituting the new variables in 
equations (^) and (^), we obtain 

K = + sgn{e,)- il - |f j (1 - C^g) + — (CgSn + f + 0{e,,e,) . (15) 

where sgn{ec) is it, for Cc positive/negative. Whereas the scaling is useful in restricting dependence 
to one parameter, B, the dynamics still concerns two coupled degrees of freedom, and no further 
reduction is, in general, possible. However, K is algebraically simpler to handle than H and we 
use it to understand the results of numerical integrations of the full equations of motion (eqns. 12). 
The dynamics of these highly eccentric orbits differs in the oblate and prolate triaxial potentials. 



4-2.1. Oblate triaxial potential 

Figures 5a and 5b are Poincare sections, strobed in g and h respectively; it is important to 
note that the iz — h section is taken at g = it /2 (mod27r), but the precise value of the strobing 
in h, for the i — g section is immaterial. Figure 5a is very similar to Figure 3a, with the notable 
addition of elliptic fixed points at /i = and vr, located at non zero values of \iz\- Moreover, 
Figure 5b shows them to be triaxial versions of the axisymmetric fixed points, discussed in § 3 (see 
Figure 2a). Recall that, in the axisymmetric case, the fixed point orbit corresponded to the steady 
rotation, about the z-axis, of an inclined, rigid Keplerian ellipse. The main effect of triaxiality is 
to introduce oscillations into quantities that were steady in the axisymmetric case. A question of 
interest is whether these oscillations change the shape of the orbit in a manner that reinforces the 
shape of the potential. Figure 6 shows the behaviour in time of g, h, i and e, where it may be seen 
that the librations of g (about 7r/2, or 37r/2) occur at twice the frequency at which h precesses. 
Physically, this can be thought of as an m = 2 perturbation of an axisymmetric problem. Below 
we provide an analysis of this perturbation, and derive expressions the quantitative effects of 
triaxiality on various quantities of interest, such as the energy at which the FPO appears, the 
amplitudes of oscillation induced in i, 1^, and g. 

Let us rewrite K (eqn. [T5| ) as the sum of two terms, Ki and K2; Ki contains terms that are 
independent of h, and K2 consists of the remaining h dependent terms. 

K = K1+K2, 
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2 4 





/2 / 



= —( —C'gC'2h. + ^5^^2/1 + y5'2g5'2/i 1 . (16) 



Ki may be viewed as an effective axisymmetric system that reduces to the actual axisymmetric 
case as 5^0. Clearly, I 2 is an integral of motion for this system, which we henceforth denote as 
Izf^. The Hamiltonian Ki admits elliptic fixed points at 



50 = ^ , -y , 



3 V 2 



1/4 

(17) 



for which the nodal angle h circulates uniformly in the retrograde sense at a frequency, 
ho = —-^^15(1 — B /2), which is independent of the value of Izq- When B equals zero, this 
expression coincides with the /iq of eqn (P), upto 0(ec). The fixed point value of the energy, i^ig 
is fixed by the value of the parameter Iz^. The minimum value of energy at which such a fixed 
point is admitted is then obtained by sustituting the above fixed point values of / and g in the 
Hamiltonian Ki, and also noting that the maximum value Iz can take is equal to /. Thus, 

Ki^.r. = \{B-l): (18) 
from which we obtain (using equation (0)), 

i?o = ^ + |ec|i^i_ = \{l-e, + eb). (19) 

The Hamiltonian K2 acts as a perturbation on this effective axisymmetric system. Below we 
present a first order analysis of the perturbed fixed point orbit. To this end, let {li,gi, Iz^, hi) be 
the perturbations to the fixed point values of the effective axisymmetric system, {Iq, goJzQ, ho). 
Recalling that {lo,go,lzo) are constant in time, and ho = — •\/15(l — B/2)s, we note that the 
observed oscillations in the values of the dynamical variables (see Figure 6) will be described by 
(hJzngiihi). Hamilton's equations of motion, given by. 



ds \dldg )^ \dlzdg j ^ \ dg^ )^ \ dg 



'^ds - [ dl'^)^ + [dlzdl ),^^\ dgdl J 

dlzj_ _ 
ds 



dhi 
ds 








(20) 
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admit the following solutions, 




/ cos2wftS, 



sin 200 hS 



sin 2uJhS , 



(21) 



where, Uh = — \/15(l — B/2). Thus it is clear that the frequency of the small amplitude 
oscillations is twice that of the mean rotation of h. Moreover, I and Iz oscillate in phase, whereas 
the oscillations of g about its fixed point value and those of h about its mean rotation are in 
phase. Note that the amplitude of these oscillations decreases to zero as B — 0, as expected for 
the axisymmetric case. Table 1 compares the predictions of our first-order theory with the results 
of numerical integrations of the exact equations of motion, for two different values of B. The 
smaller the value of B, the better is the match between the observed and predicted values. 

Figure 7 shows the projections of one such FPO along the three symmetry planes defined by, 
the long and the intermediate axes (X-Y plane), the intermediate and the short axes (Y-Z plane), 
and the long and the short axes (X-Z plane). The orbit can be viewed as a highly elongated 
ellipse, rotating about the minor axis (here Z-axis) in a nonuniform manner. The orbit makes a 
nonzero inclination that exhibits a small libration about its mean value. Corresponding to these 
librations, is a deformation of the ellipse itself, resulting in periodic changes of its eccentricity. It 
is apparent from the figures that the orbit fills a conical region in space. The orbital extent is 
maximum along the long axis, and minimum along the short axis. 

The orbit corresponding to the separatrix of Figure 5b, is shown in projection in Figure 8. 
This is seen to fill a rounder region of space than the fixed point orbit. It is interesting to note 
that both orbits (shown in Figure 7 and Figure 8) have shapes that appear to largely reinforce the 
triaxiality, as well as the oblateness of the potential. 



Figures 9a and 9b show Poincare sections, strobed at h = and g = 7:/2 respectively for a 
triaxial potential, whose prolateness is small (c<^ larger than, but close to unity). Both sections 
show FPOs surrounded by islands of stability. Recall that prolate axisymmetric potentials, with 
1 < c<^ < 2 do not admit elliptic fixed points, for any value of the energy (compare Figure 9a 
with Figure 2e). Hence this fixed point is a purely triaxial phenomenon, which is conveniently 



4-2.2. Prolate triaxial potential 
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Table 1: The observed and the predicted values of £"0, and the relative amplitudes of the fluctuations 
in £, iz and g for the fixed point orbit in oblate triaxial potential. Note that smaller the value of 
B, better is the match in values. 

^ Eq {h/h)max Q'Zi/ho)max (ffl/go)maa 

0.118 0.127 

Observed 2.334 0.548 0.031 0.059 

Predicted 2.312 0.541 0.031 0.052 



0.239 0.254 

Observed 2.356 0.170 0.065 0.120 

Predicted 2.331 0.110 0.067 0.110 



described in the £z — h surface of section. Now compare Figure 9b with Figure 3b. Both the 
configurations have the same triaxiality, but different energies. Note that the stable fixed point at 
^2 = in Figure 3b has changed its character to an unstable fixed point, as seen in Figure 9b. The 
latter figure also shows emergence of two stable fixed points, which were absent at low energies. 
It should be emphasised that these fixed points appear only above a certain minimum energy, Eq. 
The occurance of these fixed points can be explained by adhereing to an analogy with parametric 
resonance. 

Consider the Hamiltonian K, given in equation (]T5|)p|. Expanding K about the fixed point 
{lz,h) = (0,0), and keeping terms upto quadratic order in and /i, we write 

K = Kg + Kh, (22) 
= -g/^ + ^(1-^2,)), 

Kg and Kh may be viewed as the two coupled oscillators. Kh is quadratic in the small quantities, 



^ This choice of the Hamiltonian aheady implies that one is looking at potentials with small \ec\ 
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Iz and h, hence\Kh\ <C \Kg\, and we may make the approximation, 

dg ' dl ■ ^^"^^ 

These equations are identical to those of a pendulum, and are readily solved. Of the two possible 
kinds of motion, libration and circulation, it is the latter that turns out to be relevant for our 
problem. We now regard I and g as given functions of time in K^, and solve the following 
equations, 

to obtain lz{t) and h{t). Thus the dynamics has been reduced to a master-slave system, where 
the iz ~ h oscillations are driven by the i — g motion. As the energy of the system is increased, 
the forcing frquency of the Kg pendulum changes, leading to a parametric instability in the iz — h 
system. This renders the iz — h fixed point hyperbolic, simultaneously giving birth to two new 
stable fixed points, as shown in Figure 9b. The separation between the new, stable, fixed points 
increases with energy. Below, we estimate Eq by considerations appropriate to the onset of this 
parametric instability. We note that Kg is nearly constant, and eliminate / in favour of Kg and g. 



This is substituted in the expression for K^ in equation (22), to obtain 



\^ 2 {AKg + 5{l - C2g)) J V 4 y y _2(4Kg + 5il-C2g) 



Because the angle g circulates, we average the Hamiltonian, given in equation (|25|) over g, to 
obtain. 

The cross term in Iz — h has dropped out, leaving behind a simple harmonic oscillator, whose 
frequency of small oscillations is. 



'^Bil + B) l-J-^]. (27) 



^ 2 ^ V V^ + 2^V 

This should be compared with the frequency of circulation of the Kg pendulum 



-2A[Kg + ^]. (28) 



We expect parametric resonance when ujg ~ 2uJh. Imposing this condition on equations (27) and 
(|2^ , gives us an fifth order, algebraic equation for Kg, whose one real root gives us a prediction for 
the minimum energy, Eq^^ = (5/2) + led-ftTg, at which the {iz = 0, /i = vr) fixed point goes unstable, 
and bifurcates into two new elliptic fixed points. Table 2, given below, compares Eq^^, with the 
minimum energy, Eq, obtained from numerical integrations of the full equations of motion. 



Figure 10 shows projections of one such FPO on the principal planes. The orbit reinforces the 
prolateness, as well as the triaxiality of the potential. 



Table 2: The observed and the predicted values of the minimum energy Eq for the fixed point orbit 
for various prolate triaxial potentials. The superscripts obs and pre refer to the observed and the 
predicted values. 



B 



T 




0.523 



0.670 



2.36 



2.39 



0.354 



0.754 



2.31 



2.34 



0.269 



0.804 



2.27 



2.30 



4-2.3. Triaxial potential, with large prolateness 



When Cip > 2, the prolate, axisymmetric potential admits a fixed point orbit, which survives 
the introduction of some triaxiality. Figure 11 is a Poincare section strobed at /i = 7r/2 for a 
prolate triaxial configuration with c^p = 2.5 and = 0.99. This figure is very similar to Figure 2b. 
The projections onto the principal planes reveals that the orbit strongly reinforces both the 
triaxiality and the prolateness of the potential. Figure 12 shows real space projections of such a 



In the nuclear regions of galaxies, stars move in the combined gravitational fields of the 
central, supermassive black hole (BH), and the self-gravity of the cluster of stars. When the 
cluster is triaxial, the only integral of motion that orbits respect is the energy, and orbital structure 
is notoriously difficult to understand. We have applied the orbit-averaging method of ST99 to 
study stellar orbits that lie within the radius of influence of BH, where the gravitational potential 
of the cluster is a small perturbation to the Keplerian potential of the BH. Averaging over the 
orbital phase of the fast, nearly Keplerian motion promotes an additional integral of motion: for a 
time independent perturbation, this integral is the semi-major axis. Hence the problem reduces 
in dimension to the dynamics of the two precessional degrees of freedom. For this first effort at 
addressing the effect of triaxiality, we modeled the potential of the cluster by a triaxial, harmonic 
potential (the averaging method itself applies to general potentials); the rationale for this choice 
is discussed in some detail in the Introduction. We formulated the equations of motion using 
Delaunay variables (which are the action-angle variables of the Kepler problem), and integrated 
orbits numerically; appropriate Poincare sections provide global information on the dynamics. 



FPO. 



5. 



Discussion 
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The averaged dynamics of the axisymmetric case is completely integrable, and the orbital 
structure can be easily analysed. Both oblate and prolate axisymmetric potentials admit resonant 
orbits that parent families of their own. See §3 for a detailed listing of the properties of these 
orbits. Kozai (1962) discussed a similar problem in the context of the motion of asteroids 
perturbed by Jupiter. More recently, Holman, Touma and Tremaine (1997) have also considered 
the axisymmetric problem in the context of the dynamics of a planet orbiting the star, 16 Cyg B, 
being perturbed by the distant stellar companion, 16 Cyg A.^ 

In the triaxial case, the averaged dynamics admits two integrals of motions; the energy, 
which is exactly conserved, and the semi-major axis, which is secularly conserved. We undertook 
a systematic study of the nearly spherically symmetric case which, however, was quite highly 
triaxial. Chaos appeared to be, suppressed, and we identified appropriate effective third integrals 
in the several cases we studied. This could be a result of the reduction in degrees offered by 
averaging over the fast Keplerian motion. Averaging will be an increasingly invalid procedure as 
the size of the orbit approaches r^; it will not surprise us if the three dimensional dynamics is 
significantly chaotic for these larger orbits. 

Our study of the averaged dynamics of the triaxial case has turned up several resonant orbits, 
many of which appear to reinforce the triaxiality, as well as the oblateness/prolateness of the 
potential. Below we provide a brief travel guide. 

1. Polar orbits in prolate, triaxial potential : Projections on principal planes given in Figure 4. 
Poincare section in Figure 3b. Discussed in § 4.1. Orbit fills elliptical annulus in the plane of the 
intermediate and long axes. 

2. Resonant orbit in oblate triaxial potential : Projections on principal planes given in Figure 7. 
Poincare section in Figures 5a and 5b. Discussed in § 4.2.1. Orbit fills a hollow conical region, of 
elliptical cross section, whose axis is oriented along the short axis of the potential. 

3. Separatrix orbit in oblate triaxial potential : Projections on principal planes given in Figure 8. 
Poincare section in Figures 5b. Mentioned in § 4.2.1. Orbit fills a three dimensional box-like 
region, whose dimensions are compatible with the shape of the potential. Orbit spends a lot of 
time in the equatorial plane. 



4. Resonant orbit 1 in prolate triaxial potential : 
Figure 10. Poincare section in Figures 9a and 9b. 
small prolateness. Shaped like a butterfly. 



Projections on principal planes given in 
Discussed in § 4.2.2. Orbit present only for 



5. Resonant orbit 2 in prolate triaxial potential : Projections on principal planes given in Figure 12. 
Poincare section in Figures 11. Discussed in § 4.2.3. Orbit present only for large prolateness. Fills 



^In both cases, the perturber was modeled as a circular ring, and the tidal field was truncated at quadrupolar 
order, giving rise to an axisymmetric, harmonic potential of the form given in our equation (^), with $o <Q,bip — 1, 
and = 1/ \/2 . 
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a hollow gourd-like region. Highly reinforces triaxiality. A promising candidate for constructing 
self-consistent models. 
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Fig. 2a. i~ g surface-of-section for oblate axisymmetric potential, = 0.8 and energy E = 2.48 
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Fig. 2b. £ — ^ surface-of-section for prolate axisymmetric potential, dp = 2.5 and energy E = 0.5 
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Fig. 2c. ^ — g surface-of-section for oblate axisymmetric potential, = 0.7 and energy E = 2.2 
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Fig. 2d. i — g surface-of-section for prolate axisymmetric potential, c^n = 2.5 and energy E = 0.6 
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Fig. 2e. I — g surface-of-section for prolate axisymmetric potential. Cin = 1.5 and energy E = 2.3 
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Fig. 3a. iz — h surface-of-section (strobed at g = vr/2) for oblate triaxial potential, = 0.99, 
Cip = 0.96 and energy = 1.8 




Fig. 3b. iz — h surface-of-section (strobed sX g = it/2) for prolate triaxial potential, hip = 0.99, 
Cip = 1.04 and energy E = 1.8 
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Fig. 4. Real space projections of a polar orbit in prolate triaxial potential, 
and energy E = 1.8 
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Fig. 5a. — h surface-of-section (strobed at g = vr/2) for oblate triaxial potential, = 0.99, 
= 0.96 and energy E = 2.48 
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Fig. 5b. i — g surface-of-section (strobed at h = 7r/2) for oblate triaxial potential, bip = 0.99, 
= 0.96 and energy E = 2.48 
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Fig. 6. Temporal behaviour of i,e,g,h for the i — g fixed point orbit in oblate triaxial potential 
with bip = 0.99, = 0.96, and E = 2.48. The long-dashed, solid, short-dashed, and the dotted 
curves, correspond to i,e,g,h, repsectively. 
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Fig. 7. Real space projections of the i — g FPO in oblate triaxial potential, b,^ = 0.99, = 0.96 
and energy E = 2.48 
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Fig. 8. Real space projections of the separatrix orbit in the oblate triaxial case, b,^ = 0.99, 
= 0.96 and energy E = 2.48 
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Fig. 9a. i — g surface-of-section (strobed at h = 0) for prolate triaxial potential, = 0.99, 
= 1.04 and energy E = 2.3 
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Fig. 9b. £z — h surface-of-section (strobed at g = 7r/2) for prolate triaxial potential, b^p = 0.99, 
= 1.04 and energy E = 2.3 
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Fig. 10. Real space projections of the £ — g FPO in prolate triaxial potential, = 0.99, = 1.04 
and energy E = 2.3 
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Fig. 11. £ — g surface-of-section (strobed at ^ = tt/2) for prolate triaxial potential, bip = 0.99, 
Cip = 2.5 and energy E = 0.5 
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Fig. 12. Real space projections of the I — g FPO in prolate triaxial potential, b^p = 0.99, = 2.5 
and energy £^ = 0.5 



